This notebook deals with convergence testing, which is a common task and an important step in every numerical simulation technique. In the FEM-method, one typically needs to make sure that the mesh is accurate enough to represent to system for the target wavelength. Another significant parameter is the polynomial degree, or FEM-degree.
In JCMsuite, there are myriads of possibilities to controle the mesh. You can learn more about these techniques in the GeoTutorial of JCMsuite. For the FEM-degree, JCMsuite provides an adaptive algorithm to automatically chose different FEM-degrees on different elements of the mesh. You can specify minimum and maximum FEM-degrees for this procedure. The target precision, which controls the choice of the FEM-degrees, can be set using the PrecisionFieldEnergy parameter in the Refinement section of the project.
In pypmj, a class called ConvergenceTest is provided, which aims on automating the task of a convergence testing in the most general way. We will use the mie2D_extended project to showcase the approach. The meshing will be controlled using side length constraints of the involved domains.
In [ ]:
%%javascript
require(['base/js/utils'],
function(utils) {
utils.load_extensions('IPython-notebook-extensions-3.x/usability/comment-uncomment');
utils.load_extensions('IPython-notebook-extensions-3.x/usability/dragdrop/main');
});
In [ ]:
%matplotlib inline
import pypmj as jpy
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
try:
import seaborn as sns
sns.set_context('notebook')
except ImportError:
print 'You do not seem to have `seaborn` yet. You may check it out!'
For plotting, we load the prefered colormap from the configuration.
In [ ]:
cmap = jpy._config.get('Preferences', 'colormap')
We start by creating a JCMProject-instance. Make sure you configured the project catalog well, which is configured in the section Data under key projects in the configuration file, or specify the absolute path to the mie2D_extend project here.
In [ ]:
project = jpy.JCMProject('scattering/mie/mie2D_extended')
We now specify the keys for our convergence test. In contrast two a normal SimulationSet, we need two sets of keys, one for the so-called reference simulation and one for the test simulations. The reference simulation keys must induce a single simulation and its results will be used to calculate relative deviations, so we treat them as "analytical" results. Consequently, the accuracy of the reference simulation should be much higher than any of the test simulations. The ConvergenceTest-class initializes two simulation sets behind the scenes, one for the reference simulation and one for the test simulations. They will be placed in subfolders of the storage folder.
In [ ]:
keys_test = {'constants' :{},
'parameters': {'fem_degree_max': np.arange(2,6),
'precision_field_energy':np.array([1.e-2,1.e-4])},
'geometry': {'radius':0.3,
'slc_domain': np.linspace(0.1,0.4,10),
'slc_circle': np.linspace(0.05,0.3,10),
'refine_all_circle':4}}
keys_ref = {'constants' :{},
'parameters': {'fem_degree_max': 6,
'precision_field_energy': 1.e-9},
'geometry': {'radius': 0.3,
'slc_domain': 0.1,
'slc_circle': 0.05,
'refine_all_circle': 6}}
Now, the ConvergenceTest can be initialized. We will store the results in a subfolder of the current working directory for now, ignoring the storage base set in the configuration file.
In [ ]:
ctest = jpy.ConvergenceTest(project, keys_test, keys_ref, duplicate_path_levels=0,
storage_folder='convergence_test',
storage_base=os.path.abspath('tmp_storage_folder'))
The syntax for preparing and running the ConvergenceTest is prefereably identical to the one for a SimulationSet. We make a simulation schedule.
In [ ]:
ctest.make_simulation_schedule()
To run the convergence test, we need to define the processing function, which is in this case identical to the one that was used in the mie2D example. It simply reads the scattering cross section from the FluxIntegration post process.
In [ ]:
def read_scs(pp):
results = {} #must be a dict
results['SCS'] = pp[0]['ElectromagneticFieldEnergyFlux'][0][0].real
return results
Running a convergence test is nothing else than running two simulation sets, that's why you can simply pass the same arguments to the run-method. But there are two more arguments which can be specified. If save_run is True, it runs both simulation sets using the utility-function run_simusets_in_save_mode. More interesting is the run_ref_with_max_cores argument, which is set to 'AUTO' as default. As the reference simulation is a single simulation which typically needs massive computation power, it makes sense to use as many cores as possible for it. The default behavior is to detect the resource with the most available cores automatically and to use only this machine with a multiplicity of 1. For the number of threads, the product of the currently configured multiplicity and number of threads is used.
Let's make that clear. Let's configure the localhost as our only resource and set its multiplicity to 4 and the number of threads to 1:
In [ ]:
ctest.use_only_resources('localhost')
jpy.resources['localhost'].set_m_n(4,1)
print jpy.resources['localhost']
There are now 4 cores available and there is no other machine available, as we restricted the resources to localhost. The run-method will now use this machine with a multiplicity of 1 and 4 threads for the reference simulation, which is the most efficient choice. For the test simulations, the normal configuration as specified above will be used.
We now run the convergence test, 200 simulations at a time, while zipping the working directories.
In [ ]:
ctest.run(N=200, processing_func=read_scs, wdir_mode='zip')
Analyzing and interpreting a convergence test is a task which requires some experience and maybe sometimes also intuition. So the name of this section and of the analyze_convergence_results-method may be misleading in some way. The method just automates the most common step in this procedure: it calculates the relative deviation of the test data from the reference data for given result columns.
We only calculated the scattering cross section ('SCS'). In other cases, we could pass a list of result keys, i.e. column names, for the dev_columns-argument of the analyze_convergence_results-method. The sort_by-argument specifies the column of which the deviation data is used to sort the results in the end. If it is None, the first one is used.
In [ ]:
analysis = ctest.analyze_convergence_results('SCS')
analysis.head()
The returned data frame (which is also stored in the class attribute analyzed_data) contains a column called 'deviation_SCS' and is sorted by this column in ascending order. So the first row represents the simulation with the smallest deviation.
In general, the task is to find the simulation parameters which satisfy our requirements on accuracy, while causing the smallest computational cost, e.g. CPU time. As the computational costs are automatically recorded by pypmj, we have full access to that data and can plot, for example, the distribution of CPU time as a function of the SCS deviation. We color the scatter dots by the used FEM degree.
In [ ]:
fig = plt.figure(figsize=(14,4))
ax=plt.gca()
analysis.plot.scatter(ax=ax, x='deviation_SCS', y='AccumulatedCPUTime',
c='fem_degree_max', s=50, logx=True,
cmap=cmap)
plt.autoscale(tight=True)
# Save for later:
xmin, xmax = ax.get_xlim()
ymin, ymax = ax.get_ylim()
plt.show()
We observe that a higher maximum FEM-degree leads to a better accuracy. However, the relative deviation seems to saturate already for a fem_degree_max of 4 if other parameters are chosen well. This could be due to a smaller side length constraint (SLC). We could check this by calculating an average SLC and use this for the coloring. But we still have far to many parameters to make sense of the data by just looking at a single plot.
In a more advanced step, we can try to divide our data into subsets which have a constant maximum FEM degree and a constant precision field energy. Each of these subsets will be plotted seperately using the plotting technique above, this time coloring by the average SLC.
In [ ]:
from mpl_toolkits.axes_grid1 import ImageGrid
sns.set(style="ticks")
analysis['slc_average'] = (analysis['slc_circle'] + analysis['slc_domain'])/2.
max_fems = np.sort(np.unique(analysis['fem_degree_max']))
precisions = np.sort(np.unique(analysis['precision_field_energy']))[::-1]
nrows = len(precisions)
ncols = len(max_fems)
fig = plt.figure(figsize=(15,8))
grid = ImageGrid(fig, 111, # similar to subplot(111)
nrows_ncols = (nrows, ncols), # NxM grid of axes
axes_pad = (0.25,0.5), # pad between axes in inch.
share_all=True,
label_mode = 'L',
cbar_location = 'right',
cbar_mode = 'single',
cbar_pad=0.2,
cbar_size = '5%',
aspect=None)
vmin, vmax = (np.min(analysis['slc_average']), np.max(analysis['slc_average']))
for index, ax in enumerate(grid):
row, col = divmod(index,ncols)
prec = precisions[row]
fem = max_fems[col]
data_ = analysis[analysis['fem_degree_max'] == fem]
data_ = data_[data_['precision_field_energy'] == prec]
scat = ax.scatter(x=data_['deviation_SCS'], y=data_['AccumulatedCPUTime'],
c=data_['slc_average'], s=50, cmap=cmap, vmin=vmin,
vmax=vmax)
ax.semilogx(True)
ax.set_xlim((xmin,xmax))
ax.set_ylim((ymin,ymax))
ax.set_ylabel('AccumulatedCPUTime')
ax.set_xlabel('deviation_SCS')
ax.set_title('max FEM={}, precision={}'.format(fem, prec), y=1.03)
cax = grid.cbar_axes[0]
cbar = plt.colorbar(scat, cax=cax, label='slc_average')
plt.show()
Sure, this afforded some serious plotting and it may already be close to the limit of how many parameters we can study by this technique. However, we now see different tendencies as expected:
We also spotted candidates which give a very high accuracy but small computation times, in the third column of the first row:
In [ ]:
candidates = analysis[(analysis['fem_degree_max']==4) &
(analysis['precision_field_energy']==0.01) &
(analysis['deviation_SCS']<=1.e-4)]
The candidate with the smallest accumulated CPU time has the following parameters:
In [ ]:
relevant_cols = ['AccumulatedCPUTime', 'TotalMemory_GB', 'Unknowns',
'fem_degree_max', 'precision_field_energy', 'slc_circle',
'slc_domain', 'deviation_SCS']
optimum_1 = candidates.ix[candidates['AccumulatedCPUTime'].argmin()]
optimum_1[relevant_cols]
The seaborn package provides a lot of functionality to analyze data sets with lots of attributes, i.e. categories. This can be an enormous advantage and save much time. As you have seen above, plotting can be tedious in such cases and we even helped ourselves with the average SLC, which might not work out well in all cases.
Let's see what seaborn can do for us. We restrict our data to those simulations which have a deviation smaller than 0.1. (we only round the SLC values here for better readability in the plot below)
In [ ]:
subset = analysis[analysis['deviation_SCS'] <= 0.1].copy(deep=True)
for slc in ['slc_circle', 'slc_domain']:
subset.ix[:,slc] = np.round(subset[slc], 2)
The factorplot-method can produce grids of plots by only specifying the keys that we want to use for generating the rows and the columns. Each plot will be a swarmplot, showing us the deviations grouped by the FEM degree and coloured by the SLC of the domain. We will have a column for each value of precision_field_energy and a row for each slc_circle.
In [ ]:
sns.set(style='darkgrid', font_scale=1.2)
grid = sns.factorplot(x='fem_degree_max', y='deviation_SCS', hue='slc_domain',
col='precision_field_energy', row = 'slc_circle',
data=subset, kind='swarm',
palette=cmap, margin_titles=True, size=2, aspect=2)
plt.show()
We could also use the inverse of the product of deviation_SCS and AccumulatedCPUTime as a performance measure, as this is the quantity we want to minimize.
In [ ]:
subset['performance'] = 1./(subset['deviation_SCS']*subset['AccumulatedCPUTime'])
grid = sns.factorplot(x='fem_degree_max', y='performance', hue='slc_domain',
col='precision_field_energy', row = 'slc_circle',
data=subset, kind='swarm',
palette=cmap, margin_titles=True, size=2, aspect=2)
plt.show()
If we look for the performance maximum in this subset of our data, we find the following result:
In [ ]:
print 'The new optimum:'
optimum_2 = analysis.loc[subset['performance'].argmax()]
optimum_2[relevant_cols]
If we compare to the previous result, we find that we now found a candidate with a smaller deviation at smaller CPU time. This is achieved by a much smaller value for the SLC of the circle, which seems to have the dominating effect. The domain does not need an as fine grid, which is expected because of the smaller permittivity. When looking at the second figure that was created using seaborn, this behavior can be seen very clearly.
In [ ]:
print 'The old optimum:'
optimum_1[relevant_cols]